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ABSTUAt.T 


This  report  describes  several  computerized  multiple-variable,  multiple- 
response  optimization  procedures  developed  for  use  in  connection  with  simu- 
lation models.  The  four  optimization  procedures  evaluated  in  this  research 
were  as  follows: 

1.  Box's  complex  search  method  [8]; 

2.  A sequential  first-order  response  surface  approach  12,3] ; 

3.  A second-order  response-surface  experimental  design  followed 
by  mathematical  optimization; 

4.  A Box  complex  search  followed  by  a second-order  response 
surface  optimization  (methods  1 and  3 combined). 

Each  of  these  optimization  methods  involved  an  adaptation  of  the  original 

procedure  on  which  it  was  based,  in  order  to  accommodate  multiple  simulation 

responses  p . , j-1,  ...,  m. 

1 

These  four  optimization  procedures  were  evaluated  with  three  computer 
simulation  models: 

i.  A stochastic  inventory  model  [15,16]; 

ii.  A Lank  duel  model  [21]; 

iii.  A minefield  evaluation  model  [1]. 

Exploratory  work  was  also  performed  with  a simulation  model  of  a missile1 
attack  on  a surface  fleet  [18]. 

This  research  has  produced  workable  optimization  procedures  for  inter- 
facing with  simulation  models  of  naval  operat  ions,  but  it  is  concluded  that 
tb  is  interface  must  be  managed  on  a "customized"  basis  for  each  simulation 
model  due  to  the  variety  of  approaches  employed  in  naval  systems  modeling. 


INTKOm'Ci  ION 


The  design  ;ind  analysis  of  complex  systems,  especially  those  involving 
probabilistic  or  stochastic  elements,  often  necessitates  the  use  of  digital 
computer  simulation.  A problem-solving  procedure  for  defining  and  analyzing 
a model  of  a system,  simulation  affords  the  analyst  the  opportunity  to 
evaluate  complex  systems,  (a)  without  constructing  them  if  they  are  proposed 
sy:;U ms,  (b)  without  interrupting  their  operation  if  they  are  actual  systems, 
or  (c)  without  damaging  or  destroying  them  if  the  purpose  of  simulation  is 
to  test  the  effectiveness  of  the  system  in  a hazardous  operating  environ  nt. 

Naval  opera t ions  present  especially  challeti;  ing  problems  for  the 
simulation!.'; t . The  purpose  in  simulating  naval  operations  would  very  likely 
be  to  evaluate  the  effectiveness  of  certain  offensive  or  defensive  fleet 
configuration,  and  tactics.  Large  numbers  of  controllable  input  variables, 
uncontrollable  input  variables,  and  measures  of  system  effectiveness  are 
often  involved  in  realistic  models  of  naval  operations.  This  not  only 
complicates  model  development , but  lieu  pers  experimentation  with  the  model 
as  well. 

Simulation  can  lie  defined  as  the  establishment  of  a mathematical- 
logical  model  of  a system  and  the  experimental  manipulation  of  that  model 
on  a digital  co;  putcr.  This  paper  concentrates  on  the  latter  aspect, 
particularly  upon  means  of  determining  an  "optimal"  design  through  experi- 
mentation with  a digital  computer  simulation  model.  (The  word  "optimal" 
is  used  guardedly  here,  > inco  the  randomness  inherent  in  computer  simulation 
i trad!  ts  th<  cl  tsical  notion  of  an  optimal  solution.)  This  research 
he  explored  procedures  for  combining  computer  simulation,  response  sur- 
t thodol  Leal  pr<  ra  ling  to  seek  an  estimated  "opt  i > ' " 


I 

I 


solution.  These  technique;,  are  demonstrated  and  compared  using  simulation 
models  of  (1)  a stochastic  inventory  system  (2)  a task  duel,  and  (3)  mine- 
field evaluation. 

BACKGROUND 


A computer  simulation  model  can  he  regarded  as  a "black-box",  as 
illustrated  in  Figure  1,  in  which  values  for  ji  controllable  input  variables 
Xj , i=l,...,n  are  combined  in  some  manner  to  produce  values  for  a set  of  m 
output  or  response  variables  q . , j=l, . . . ,m.  These  responses  are  measures  of 
merit  or  effectiveness  for  the  system  being  modeled.  It  is  typically  true 
that  the  system  responses  are  affected  by  a set  of  uncontrollable  factors, 
zy,  k-1 , . . . , p . 

The  execution  of  the  computer  simulation  model,  cither  for  some  fixed 
simulated  time  t or  until  r.  realizations  of  the  j-th  syst  2m  response  have 
been  obtained,  produces  m simulated  responses. 


y.=g.  (xJ,...,Xn|„1,...,Zp)+  e.,  j=l , . . . ,m 


0) 


The  quantity  y^  is  an  estimate  of  the  true  system  response  n ^ , and  is 

found  from  the  relation 
r 

j 

y.  = ~~  l r>ni  > 3=1,..., m (2) 


t*  " -f  ^ 

j rj 


where  T, . „ is  the  individual  value  of  the  j-th  system  response  at  the  5,-ih 
realization  and  r.  is  the  number  of  realizations  of  this  response.  Equation 
(2)  suggests  that,  for  a given  execution,  or  simulation  trial,  the  r ^ , 


j=l,...,n  need  not  be  the  same.  The  simulation  trial  also  produces  an 


unbiased  estimate  of  the  variance  of  the  j-th  system  response;  that  is. 


4 


r.i 


<vf>  2i  (^~y3)J 


(3) 


or  more  conveniently 
1 


2 
s . 
J 


(r  .-1) 


r . 

3 2 ; 

> f - r . v . 

£=1  3 '3 


(3a) 


Fishman  f 1 2 ] and  Kleijnon  | 19]  describe  other  procedures  lor  computing  s',  in  the 

case  of  serial  correlation  of  the  relationships  C . . 

J*- 

The  structure  of  equation  (1)  must  he  carefully  examined  in  order  to 
develop  the  basic  framework  for  the  procedures  presented  here.  It  is  stated 
that  y.  is  an  estimate  of  the  true  response  n.  ; hence, 


y . = n . + c . 

3 3 .1 


(A) 


The  Ej  term  reflects  the  random  variation  inherent  in  the  probabilistic 
system  under  study,  and  can  he  regarded  as  sampling  error.  The  uncontroll- 
able factors,  the  z, , k=l,...,p  terms  in  equation  (1),  can  be  viewed  as 

K. 

contributing  to  the  random  variation  e..  Tliuy  equation  (1)  can  be  simplified 
to 


y.  = g (Xj  , . . . ,xn)  + Cy  j = ],..., m (5) 

Tlic  estimate  y is  therefore  a random  variable  with  mean 

j 

K(y  j)  = nj  =g.(xi,...,xn),  j=]  , . . . ,m  (6) 

and  variance 

Var  (y , ) =Var(r.^,  j = l m (7) 

This  statement  suggests  that  the  true  response  n . is  a point  on  an  unknown 
(ii+) ) -dimensional,  hypei  :.ur  face  g (X)  at  a given  point  X -•  (x^ , . . . .x^)  and 
that  the  simulation-generated  response  v.  is  a variate  from  a probability 


r 


distribution  about  the  true  value-  11 . at  this  point  X. 

The  response  surface  methods  employed  in  this  methodology  require  the 

assumption  and  the  random  error  c.  lie  normally  distributed  with  mean  zero 

, • 2 
and  variance  that  is. 


E(y . | x , . . . ,x  ) = n . 
j 1 11  j 


j=l. • • • ,rn 


Var(y. Ixj, . . . ,xn)  = o , j=l, 


,m 


(8) 

(9) 


It  is  also  required  that  Yar (y . j . ,x^)  be  the  same  at  any  point  on  the 
( n+1 ) -d im en  s i o n a 1 hypersurface. 

Note  that  there  are  m distinct  (n+1) -dimensional  hypersurfac.es,  each 
having  its  own  characteristic  random  variation.  The  problem  at  hand  is  to 
have  tiic  computer  simulation  model  which  produces  these  m responses  be  either 
automatically  or  interactivel;  controlled  by  an  "optimization"  procedure, 
as  illustrated  in  figure  2.  In  general,  the  "optimization"  procedure  must 
combine  the  experimental  features  of  response  surface  methodology  with  the 
logical/numerical  procedures  of  mathematical  programming. 

PROBLEM  FORMULATIONS 

We  shall  consider  two  basic  approaches  to  formulating  the  problem  of 
optimizing  a multiple-variable,  multiple-response  simulation  model.  One 
approach  is  the  familiar  constrained  optimization  formulation  in  which  one 
of  the  simulation  responses,  say  n , is  to  be  maximized  cr  minimized, 
subject  to  maintaining  the  remaining  m-1  responses  within  prescribed  bounds. 
The  second  approach  is  the  multiple-objective  formulation,  in  which  the  m 
responses  are  either  weighted  to  form  a single  objective  or  treated  in  a 
manner  akin  to  goal  programming.  Each  of  these  two  formulations  is 
described  in  the  following  section:;. 


1011 


Cons t r.'i  i si  .1  Opt  in'  i / 


Under  the  coma  rained  optimisation  approach , the  problem  is  stated  as 


Max  (or  Min)  r\  = g (x  ) 

1 1 1 n 


subject  to  the  constraints 


a . < x . < c . , i=l , . . . ,n 
x — l — i 


"j  ’ Sj(!:i 


dj  ’ » 


The  constraints  expressed  in  (11)  are  bounds  on  the  controllable  input 


vari  1 x and  are  typically  ! j Lori.  The  bounds  (11) 


genera]  y form  th<  ki  ent  il  rc  ;ion  prior  to  conducting  si 


trials.  In  contrast  to  that,  the  response  functions  g (Xj .... .x^)  in  (1?) 


are  not  usually  known  a pjmlor i and  hence  the  responses  n . must  bo  esti- 


mated experimentally  via  simulation.  Thus,  simulation  trials  performed  t 


points  satisfying  (11)  may  yield  responses  violating  (12) 


To  complicate  i. utters  even  more,  the  random  error  c.  can  lead  to 


erroneous  decisions  relative  to  feasibility  with  respect  to  the  constraint  - 


in  (12).  The  same  is  true  relative  to  the  objective  function  in  (10), 


That  Is,  one  simulation  trial  can  appear  to  represent  an  improvement  over 


another  when  the  true  response  at  this  particular  set  of  values  x ^ x 


does  not.  These  difficulties  can  be  countered  through  variance  reduction 


techniques,  as  discussed  by  Fishman  [12]  and  Kleijncn  [19] 


Multipl  ] ■ | ctl*  0 ' : is  t i.oi 


<>  approach  to  a multiple-obj ectiv<  f<  rmul  ition  is  to  assi  n» 


, ' n to  tl  r<  1 and  fo  ingli  tion 


f 


I 


The  bounds  (.1.1)  still  apply,  so  that  the  problem  remains  one  of  constrained 

optimisation,  but  one.  in  which  the  entire*  feasible  region  is  known  a priori 

The  weight  w , , j-l , . . . ,r.i  are  typically  assignee!  through  the  subjective 

judgment  of  tin-  decision-maker  in  the  system  being  simulated.  These  weight 

are  usually  normalized,  so  that 
m 

l w =1  (14) 

j=l  ' 


One  frequently  encounters  the  situation  in  which  certain  of  the  responses 
rp.  are  to  lie  : .ximized  and  otlior  minimized.  This  case  is  handled  by  maxi- 
mizing the  negative  of  those  functions  which  are  to  be  minimized,  so  that 
tlie  objective  function  in  (13)  is  rearranged  to  the  form 

s m 


Max  VI  ■-  } w.g.  (x,  , . 
j-l  1 J 1 


••.x,)  - l 


j-s+1 


vy..  xn) 


(15) 


where  s functions  are  maximized  and  m-s  functions  are  mini  lined. 

A second  approach  to  Lhe  multiple-objective  formulation  is  one  which 
casts  the  objcf  Live  function  as  a "utility  function";  that  is. 


subject  to  the  bounds  in  (11).  The  formulation 
of  that  in  (lb),  in  which  Ulgj  (Xj  , . . . .x^)  ] is  a 
Montgomery  and  bettenc.ourt  [31]  discuss  various 
multiple  objective  optimization  problem,  as  well 
it;  solution,  and  demonstrate  its  application  to 
s 1 mu  1 at  ion. 


06) 

in  (15)  is  a special  care 
linear  additive  function, 
formulations  of  the 
as  several  approaches  to 
mill  tipi  e-response 


Another  multiple-objective  optimization  formulation  is  that  called 


goa  1 programm i ng . This  procedure  is  ini  t ialod  hy  cstahl  i sh  i ng  a set  of 
goa  1 s in  terms  of  the  m system  responses.  These  goals  are  expressed  as 


C.=S.(V...,Xn),  j=l , . . . ,m 

Each  goal  must  have  an  associated  right-side  value  d.;  that  is, 


(17) 


GJ  = f’’j  (nj V 


d.,  j=l 


(18) 


With  a slight  modification,  each  goal  can  he  expressed  as  an  equality 


C.  =g.(>:r...,xn)  + «.  -p.  =d.J<L,...,m 


(19) 


where  n.  is  a negative  deviation  from  d . , and  p.  is  a positive  deviation. 
J J J 

Either  n.  or  p.  must  he  zero  in  any  given  solution,  and  both  could  he 
J .1 

zero.  Next,  each  of  the  n goals  C.  is  assigned  to  a priority  level  Pj  , 
k=l,...,t,  where  F ^ represents  the  highest  priority  and  Pf  the  lowest. 

For  any  goal  falling  within  a given  priority  level  Pj  , the  decision-maker 
eigl  oai  re ! at  ive  to  another.  The  i tnal  iti  p Ln  problem  formula- 

tion is  to  combine  these  several  levels  of  goals  into  an  achievement 
function  which  has  the  form 


{'’l^V^P’  r2('n2’p2^ 


(20) 


This  achievement  function  is  simply  an  ordered  C-vector.  Its  structure  is 
predicated  on  one  of  the  following  procedures  for  achieving  the  j-th  goal: 
(a)  To  equal  nr  exceed  d.,  minimize 
(h)  To  equal  or  he  less  than  d.,  minimize  p. 

(c)  To  equal  d.,  minimize  (n . + p.) 

A solution  (y’‘,...,x*)  is  con.  idered  optimal  if  the  corresponding  value  A,v 


is  the  same  as  or  preferred  to  any  other  value  A.  Therefore,  the  general 


JO 


goal  programming  problem  i..  to  find  no  ns  to  minimize  the  ord<  red 


vector  (20)  such  that  tin' 

goals  (19)  are  sat  isfivd 

and 

X . > 

1 — 

0,  1=1,. 

. . ,n 

n . > 

J - 

0,  i=l,. 

. . , m 

(21) 

pj  - 

0,  j=i,. 

. . ,ra 

The  functions  g 

( V 

J 1’  ’ 

< ) in  (19)  are  generally 

unknown , 

but  arc  usually 

assumed  to  be  nonlinear.  Any  technique  proposed  for  solving  this  problem 
in  the  simulation  domain  must  provide  experimental  estimates  of  these  unknown 
functions,  as  well  as  a mathematical  procedure  for  optimization.  Moreover, 
the  experimental  observations  are  produced  via  simulation  - etch  simulation 
trial  at  a point  Xj,...,x^  produces  m responses  v.=l,...,m.  Biles  [4] 
has  described  the  application  of  nonlinear  goal  programming  to  the  multiple- 
response  simulation  problem,  based  on  techniques  proposed  by  Tgnizio  [17]. 

Ol'TTMl  /ATI ON  TfCir MOTTS 

Various  procedures  have  been  applied  in  combining  optimization  and  simu- 
lation to  seek  the  "optimum"  solution  to  systems  possessing  a single  response 
n . The  multiple-response  problem  described  here  is  complicated  by  the 
necessity  to  observe  several  responses  at  once,  and  to  incorporate  these 
values  into  the  optimization  technique.  But  many  of  the  same  techniques 
that  have  been  applied  successfully  to  the  single-response  problem  can,  with 
appropriate  modifications,  be  extended  to  accommodate  multiple  responses. 
Moreover,  Lhasa  modified  procedures  arc  often  applicable  to  more  than  one  of 
the  aforementioned  formulations  of  the  multiple-response  problem. 

The  optimization  pro  cdures  described  below  fall  into  three  categories: 
(1)  direct  search  techniques,  (?.)  first-order  response  surface  methods,  and 


1 1 


( ) )i  ,i  : r fa  * pr<  . Althougl  n terow  techniques 

will  : • . ' , . . few  br  . 1 • stat  d procedures  will  be  outlined  here. 

It  should  be  rt.'.'K  ■ Iv  red  lli.it,  a 1 though  we  -av  refer  to  "optimization" 
lechniqui  the  ela-  .ical  notion  of  an  "optimum"  solution  is  inapplicable  duo 
to  the  pro  once  of  the  sampling  error  t.  associated  with  each  response  vari- 
able n. • Rather  we  shall  seek  a solution  which  hopefully  lies  close  to  the 
true  so] ut ion.  In  a rore  fornal  sense,  we  might  state  that  we  are  to  some 
d,  :;re.  confident,  say  90%,  that  the  true  solution  lies  within  some  interval 
about  cur  estimated  solution. 

Dire c t S'  . irc1  s Me t he d s 

Direct  search  methods  are  those  which,  applied  in  a purely  computational 
manner,  do  not  require  the  use  of  derivatives.  These  methods  progress 
through  a sequence  of  points  according  to  c,;ic  algorithm.  Typical  of  this 
class  of  optimization  techniques  are  the  pattern  search  algorithm  by  Hooke 
and  Jeevi  • [14),  sequential  simplex  search  by  Spcndley,  llext  and  lliusworth 
[26],  and  the  so-called  "complex"  search  method  by  M.  .1.  Box  [8].  In  general, 
these  direct  search  procedures  make  rapid  early  progrt  s toward  an  "optimum", 
hut  iterate  laboriously  as  a solution  is  neared.  This  is  particularly  true 
in  the  presence  of  random  error,  as  encountered  in  simulation. 

Among  the  direct  search  techniques,  Box's  "complex"  method  18)  has  been 
found  to  be  the  easily  adapted  to  a multiple-response  environment.  It  ale 
performs  better  than  any  of  the  other  direct  search  techniques  in  the  face  of 
random  crroi  and  c<  raints.  In  fact,  "<  plex  earch  is  not  at  al i 
complex,  but  derives  its  name  from  a contraction  of  the  words  "constrained 
i lex":  i i evo Ived  fro  th«  [ucntial  Imp] ex  met h >d  I 26]  and  th« 

necessity  to  deal  with  constraints.  Tiiis  report  describes  a modification 
of  Box's  method  whicl  It  especiall)  ' b 1 < I < the  multi  resp 


! 


12 


simulation  problem.  The  following  procedure  doscribi  s a generalized 
"complex"  proeodure  a;  ii  might  1h*  applied  ti>  tin  multiple-response  simulation 
problem: 

1 N 

1.  Randomly  generate  a set  of  N n + 2 search  points  X ,...,X 
satisfying  the  known  bounds  (11). 

2.  Per  form  a simulation  trial  at  each  of  these  I.'  search  points 
and  record  the  : X estimated  responses  y . , j = l,.,.,m, 

1-1,..., N. 

3.  Where  a given  point  x'‘  is  observed  to  violate  one  or  more 
constraints,  if  such  constraints  apply  with  the  particular 

problem  formulation  being  employed,  generate  a replacement 

k » k' 

search  point  X , perform  a simulation  trial  at  X , and 

k ' 

record  the  n estimated  response  at  X'  . 

4.  After  li  feasible  search  points  ha  . e been  established,  evaluate 

the  objective  function  for  each  of  these  X points.  This 

"objective  function"  might  bo  r,j  in  (10),  W in  (13),  U in 

(16),  or  A in  (20).  Among  these  N search  points,  find  the 
w 

worst  point  X';  that  is,  the  search  point  giving  the  least 

desirable  value  of  the  objective  function.  Define  X1  as  the 

w 

centroid  of  the  N-l  remaining  points.  Project  from  X 
through  X1  to  the  image  point  Xv  . If  the  known  bounds  (11) 
are  violated  by  this  move,  shorten  1 he  step  to  Xw  until  no 
violation  occurs.  Perform  a simulation  trial  at  Xw  . 

5.  Repeat  steps  3 and  4 until  a solution  (X  ,Y  ) is  obtained  which 
represents  the  best  solution  that  can  be  achieved  within  the 
available  computer  lime.  Figure  3 illustrate  Don's  complex 
search  as  applied  to  a constrained  problem.  Figure  4 gives  a 


Fipu  r /) . A multiple-response  complex  search  procedure 
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a flow  chart  for  a s iimil  at  iou/opt  in  i za  t ion  method  based  on 
the  procedure  oaf  lined  above.  Appendix  A contains  the 
program  listing  for  a constrained  optimization  formulation 
of  the  Hox  eon. <lex  search  method. 

A significant  advantage  of  complex  search  is  that,  once  an  initial 
"complex"  of  N feasible  simulation  trials  (and  perhaps  several  infeasible 
trials)  have  been  performed,  trials  are  conducted  one  at  a time  thereafter. 

The  search  ean  be  continued  as  long  as  improved  solutions  are  obtained.  If 
several  successive  simulations  are  performed  at  scattered  points  around  the 
known  experimental  re  ion  without  achieving  an  improved  solution,  however, 
the  search  can  be  terminated  and  the  best  available  solution  adopted, 

First -Order  Kespons  Surface  Methods 

First-ord.  r response  surface  methods  attempt  to  accomplish  experimentally 
what  the  "method  of  steepest  ascent"  accomplishes  computationally.  From  a 
current  point  Xk,  a designed  experiment  is  conducted  (with  a simulation 
trial  at  each  d-  ign  point)  to  estimate  the  gradient  direction  Vg(Xk). 
Simulation  trials  are  then  conducted  at  points  along  this  direction  to  a new 
point  Xk+I  which  represents  the  best  solution  obtained  along  the  direction 
Vg(Xk).  This  process  is  an  experimental  approximation  of 

Xk+1  = Xk  + XkFVg(Xk)]  (22) 

The  step  length  \k  can  be  estimated  by  a line  search  or  by  a regression 
procedure  as  dr  < 1 ibed  by  Piles  [2,31. 

Tlie  gradient  direction  Vg(X  ) is  estimated  by  placing  an  appropriate 
first-order  experimental  design,  such  as  a 2n  factorial,  2*'  ' fractional 
factorial,  or  n-dinun  tonal  simplex  design  (see  Myers  [22])  around  the. 
current  i ini  Xk.  A simulation  trial  Is  performed  at  each  point  in  the 


10 


selected  experimental  design.  From  these  N observations  the  multiple  linear 


regression  model 


y 


b 

o 


n 


+ z 


(23) 


can  be  estimated  (see  Draper  and  Smith  101).  Since  the  gradient  direction 
Vg(Xk)  is  mathematically  defined  as  the  n-vector  of  first  partial  derivatives 
of  g(X)  evaluated  at  Xk,  it  is  clear  that  Vg(Xk)  is  simply  the  n-vector  of 
regression  coefficients  exclusive  of  the  b0  term;  that  is. 


Vg(Xk)  = (b1,...,bi)'  (2/,) 

In  the  multiple-response  simulation  problem,  a simulation  trial  is  con- 
ducted at  each  design  point  in  the  selected  first-order  design  and  the 

l 

m observations  y j , j=l,...,m  are  recorded  at  each  design  point.  Multiple 
linear  regression  is  applied  separately  to  each  set  of  observations  (assuming 
independence  among  the  m responses),  producing  the  m models 

n 


y . = b . „ 
3 J,0 


+ z 


b. 

3 i 


l*i> 


j=l m 


(25) 


and  hence  the  m gradient  vectors 


Vg.(X')  <=  (b  ,...,b  )'  , j=l,  ...,m 

J J y -*■  J y 


(26) 


These  estimates  can  then  be  employed  in  any  one  of  several  optimization 
schemes  to  produce  an  improved  solution  y}' 1 ^ . A generalized  procedure  for 
accomplishing  this  improved  solution,  and  an  estimated  "optimum",  will  be 
described  later.  But:  first  it  is  necessary  to  give  attention  to  the  experi- 
mental designs  employed  to  estimate  the  gradient  vectors  Vg^ (Xk) , j=l,...,m. 

In  selecting  a first-order  response  surface  design,  it  is  usually 
desirable  to  minimize  the  variances  of  the  regression  coefficients  , i“l , . . . ,n. 
To  acre:  plisli  this  the  first-order  experimental  design  should  he  orthogonal. 


17 


That  is,  the  placement  of  the  N experimental  poinLs  (in  our  case,  simulation 
trials)  is  described  by  the  N by  n design  matrix  D,  where 


11 

• • • 

Xnl 

12 

x22 

* * * 

Xn2 

IN 

X2N 

XnN 

(27) 


Then  an  N by  n+1  matrix 

X 

is  constructed  by  placing  a 

unit 

vector  to  the 

left  of  D.  Thus, 

1 

X11 

X21  "•  Xnl 

X = 

1 

X12 

x22  - ' Xn2 

(28) 

1 

X1N 

X2N  ’ ’ 1 xnM 

It  is  usually  convenient 

to 

code  the  design  levels, 

, SO 

that 

the  following 

conditions  are  achi 

eved : 

N 

l X 

U=1 

2 

iu 

N 

i=l , . . . ,n 

(29) 

V x.  =0 
).  1U 
u=l 

If  the  actual  value  of  the  u-th  level  of  the  i-tli  variable  is  z , then  the 
corresponding  coded  values  is 


‘iu 


z , - z . 

_IU HI 

s . 

i 


(30) 


I 


L 


I 

I 

I 
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where 


\m  ( l 

U=1 


(31) 


and 


N 


si-  l, 

u=l 


(32) 


Then 


X'X 


N 0 0 ...  0 
0 N 0 ...  0 


0 0 0 


0 


(33) 


Since  the  (n+1)  - vector  of  regression  coefficients  b is  estimated  by  the 
least  squares  relation 


b = (X'X)“[  X'  y 


(34) 


where  y is  the  N-vcctor  of  response  estimates  obtained  from  N simulation 
trials,  the  variance  of  the  regression  coefficients  Ik,  i=l,...,n  is  given 


bv 


Var  (Ik  ) = o /N, 


i-1, . . • ,n 


(35) 


where  o is  the  variance  of  the  error  term  i: . Since  wc  are  interested  in 


m separate  system  response  y.,  j=l,...,m,  equations  (34)  and  (35)  can  be 


generalized  to 


b . = (X'X)  1 X'  y , j=l, . . . ,m 


(36) 


Var(b..)  = o./N,  i=l,...,n 

ji  J 


j=l , . . . ,m 


(37) 


]‘J 


I 

I 

I 

I 
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Tor  an  orthogonal  first-order  design,  Liu  results  in  (33)  — (37)  hold,  giving 
a so-called  "minimum-variance"  design.  The  2n  factorial  and  2n  * fractional 
factorial  designs  are  orthogonal  and  hence  minimum  variance.  Orthogonal 
n-simplex  designs  can  bo  easily  constructed  (see  Myers  [22().  Since  n-simplex 
designs  provide  the  minimum  number  of  design  points  needed  to  estimate  the 
multiple-linear  regression  models  in  (22)  or  (25),  and  are  hence  the  most 
"economical"  of  the  first-order  response  surface  designs,  they  are  especially 
attractive  for  the  purpose  proposed  here.  Brooks  and  Mickey  [9]  concluded 
that  n-simplex  designs  offer  the  most  efficient  approach  to  estimating  the 
gradient  direction  Vg.(X).  Figure  5 illustrates  2n  factorial  and  n-simplex 
designs . 

Biles  [2,3]  has  described  a first-order  response  surface  procedure  for 
approaching  the  constrained  formulation  of  the  multiple-response  simulation 
problem.  This  procedure  Involves  performing  a first-order  design  around  a 
current  point  Xk  to  estimate  the  gradient  direction  'g(X^)  according  to 
relation  (24).  A line  search  Is  then  pei  ;(X  ) to  estimate  an 

optimal  step  X in  (22).  As  long  as  the  search  remains  interior  to  the 
region  hounded  by  the  constraints  (11)  or  (12),  the  procedure  is  basically 
the  same  as  that  proposed  by  Box  and  Wilson  [6].  If  one  or  more  constraints 
(11)  oi-  (12)  are  encountered,  however,  Biles  [2,3]  proposes  that  the  gradient 
projection  direction  be  followed  (sec  Korea  [23]).  The  procedure  for  estimat- 
ing the  gradient  projection  direction  is  ns  follows. 

Suppose  that  at  an  estimated  boundary  point  x'',  q constraints  arc 
satisfied  as  equalities.  These  can  be  cither  the  (11)  or  (12)  constraints, 
or  both.  I.ct  be  the  n x q matrix  of  first  partial  derivatives  of  these 
active  constraints.  Thus,  I!  consists  of  the  q gradient  vectors 
Vg.(Xk),  j-1 That  is, 


?q,/Dx  3g  / 3x 

1 n q n j 


Since  g (X) , j=l,...,q  denotes  Lhc  set  of  binding  constraint  functions 

(a  constraint  (11)  or  (12)  is  said  to  be  "binding"  if  it  is  satisfied  at  the 

equality),  for  the  moment  let  f (X)  represent  the  objective  function.  Then 

Vf(x')  and  Vg^(Xk),  j-l,...,q  represent  the  gradient  vectors  of  the 

objective  and  constraint  functions,  respectively,  evaluated  at  the  boundary 

__k 

point  X . 

Performing  a first-order  response  surface  experiment  about  the  boundary 
point  Xk  yields  estimates  of  the  gradient  vectors  Vf(X^)  and  Vg^(X^) , j=l, . . . ,q 
in  the  form  of  the  vectors  of  regression  coefficients.  (If  a constraint  of 
type  (l1)  is  included  in  the  set  of  binding  constraints,  the  gradient  vector 
has  the  form  (0,0, . , . , 1 , . . . ,0,' ' , where  all  elements  are  zero  except  the  i-th 
clement  which  is  one)-  Following  the  procedure  outlined  by  Rosen  [23],  the 
gradient  projection  direction  is  given  by 

Sk  = [Vf(Xk)]  - (B'^r1  IT  [Vf  (Xk) ) (39) 

A "golden  section"  line  search  is  performed  along  direction  S until  either 

(a)  a local  "optimum"  is  found,  or  (b)  other  constraints  are  encountered, 

k+1 

This  new  point  is  denoted  X . This  procedure  is  repeated  until  the  gradient 
projection  direction  Sk  is  approximately  zero.  This  point  X"  is  taken  as  a 
"constrained  optimal"  solution.  Figure  5 illustrates  the  application  of  the 
gradient  projection  procedure  to  a constrained  optimization  problem. 

Swain  [27]  has  compared  other  first-order  response  surface  techniques, 


of  Klin  nan  and  liinmelblnu  [20]  and  Zoutendijk  [28],  to 


Rosen's  gradient  pro /,  eel ion  method  [23).  lie  found  little  difference  among 
these  techn i ques  in  Lenns  of  experimental  requirements,  and  lienee  computer 
simulation  time,  but  saw  significant  variability  in  computational  require- 


ments in  order  to  use  these  algorithms.  The  /.out end  ! jk  methods  of  feasible 
directions  [28)  require  greater  computational  effort  than  the  other  procedures. 

Biii..-.  [•']  demonstrated  both  first-order  and  second-order  approach'  to  n 
nonlinear  goal  programming  formulation  of  the  multiple-response  simulation 
problem.  Th  e ipproacl  are  based  on  Ignizi  ' iptation  [17]  of  the 
method  of  Griffith  and  Steward  [13]  to  goal  programming  and,  like  the 
constrained  procedures  mentioned  previously,  combine  si mu'  it  ion,  experimental 
design  and  mathematical  programming. 

The  following  generalized  procedure  is  followed  in  employ inp,  a first- 
order  response  surface  approach  to  the  multiple-response  simulation  problem. 

The  particular  problem  formulation  and  opt i e will  govern  ( 

precise  sequence  of  steps  in  implementing  this  procedure. 

1.  Identify  the  known  experimental  region  a •_  x^s.  e , , 1=1 , . . . ,n. 

Selecting  a starting  point  within  this  region.  With  as  its 
center,  array  an  orthogonal  first-order  response  surface  design 
within  a selected  design  radius.  Place  n = n/2  >_  2 points  at  the 
design  center  X^  (coded  as  the  0 - vector). 

2.  Perform  simul  t ion  trials  at  each  of  the  N experimental  design  points 

and  record  the  responses  y.,  jl,...,m;  ,,X.  Using  multiple 

linear  regression,  fit  linear  models  of  the  form  (23)  and  (25). 

3.  Apply  the  appropriate  mathematical  programming  technique  to  locate 
the  next  center  point  in  the  search. 

A.  Repeat  steps  1-3  until  an  "optimum"  solution  is  located,  Tt  may  he 
appropriate  to  dd  design  points  to  complete  a second— order  response 
surface  design  t«  test  this  optimum  solution.  The  procedure  for 


I 

I 

I 
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ac.coupl  ishing  this  is  described  in  the  next  sect  ion, 
r i pure  6 illustrates  a sequent  ial  first-order  response  surface  method  applied 
to  a constrained  system.  1'igure  7 gives  a flow  chart  for  a simulation/ 
optimization  method  based  on  the  procedure  outlined  above.  Appendix  b gives 
a program  listing  for  this  procedure. 

Secoi  j ' J i i tespon  e Surl  e Methods 

A second-order  response  surface  approach  to  the  multiple-variable, 
multiple-response  milation  problem  consists  of  one  or  more  repetitions  ot 
a two-stage  procedure:  (a)  the  execution  of  a computer  simulation  trial 
at  each  point  in  a second-order  response  surface  experimental  design  cover- 
ing the  known  region  given  by  (11)  and  the  use  of  multiple  linear  regression 
to  fit  second-order  regression  models  to  the  resulting  data;  and  (h)  the 
application  of  a suitable  mathematical  programming  procedure  to  obtain  a 
solution  to  the  problem  formulated  in  (10)-(12),  in  (15)  together  with  (11), 
in  (16)  together  with  (11),  or  in  (19)-(21),  In  contrast  to  the  first-order 
methods,  in  which  the  optimization  procedure  was  part  and  parcel  with  the 
experimental  procedure',  these  procedures  are  distinct  and  sequential  in  the 
proposed  second-order  approaches. 

The  first  step  in  the  second-order  approach  is  to  identify  the  range  of 

each  Input  variable.  A safe  strategy  is  to  cover  the  entire  known  region 

a.  < x < c.,  i=l n with  the  first  (and  possibly  only)  experimental 

l — i — l 

design.  If  we  let  a,  denote  the  radius  of  the  n-dimonsional  hvpersphere 

i 

within  which  the  design  points  are  contained,  then 


“i  = ~ i=l,  • • . ,n 


(40) 


is  effectively  the  maxinun  radius  we  could  construct.  It  is  convenient  to 
adopt  the  coding  convent  ion  expressed  in  (29)  - (30),  hut  choosing 
in  such  a way  that  <.  satisfies  (40).  Myers  [22]  describes  lids  coding 
process. 
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Figure  7. 


Flow  Chari  of  Gradient  Projection  Search 
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Tlu>  second-order  fitted  response  surlace  has  the  form 


+ ) h.x:  + y b..x.^  + ) y 

O V . 1 1 .L  . 11  1 . 1 . 


i=l 


i=l 


i~l  j=l 


b. . x ,x  . 
ij  i ] 


(41) 


whore  y is  the?  estimate  of  the  true  response  n at  a given  value 

X = (x  x ) and  the  b.  and  b.  are  regression  coefficients  in  the  fitted 

in  l ij 

model.  Since  we  must  estimate  ri  separate  response  relationships,  equation 
(41)  is  modified  to 


yk  = bk,o  + ^=1  bk,i*i  + lml  bk,iixi  + bk,ij  xiXj 

i * j 

k = 1, . . . ,m 


(42) 


Given  the  independence  of  the  m responses,  these  m regression  equations  can 
be  estimated  independently  from  a set  of  1.  >_  (nil)  (n+2)/2  data  points  ob- 
tained by  performing  a simulation  trial  at  each  point  in  a second-order 
response  surface  design. 

An  experimental  design  employed  for  the  purpose  of  estimating  the 

regression  coefficients  in  (42)  must  contain  at  least  as  many  design  points 

as  there  are  coefficients  b.  and  h..  in  the  fitted  model,  of  which  there 

i ij 

are  (n+1)  (n+2)/2.  Because  of  the  non-linearity  of  (42),  the  experimental 
design  must  also  have  at  least  three  levels  of  each  controllable  variable 
x.  , 1 = 1 , . . . ,n.  It  is  also  desirable  to  have  a design  which  is  rotatable; 
that  is,  the  predicted  response  y at  some  point  X is  a function  only  of  the 
distance  from  the  design  'enter  to  X and  not  a function  of  the  direction. 

The  most  widely  used  design  for  fitting  a second-order  model  is  the 
central  composite  design,  shown  in  Figure  P for  n=2  and  n~ 3 . These  designs 
consist  of  a 2n  factorial  (or  suitable  fraction  thereof),  augmented  by  2ti 
axial  points  and  k center  points.  A central  composite  design  can  he  made 
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rota l able  by  proper  choLce  of  «,  the  distance  of  the  axial  points  from  the 
design  center.  With  Clio  proper  choice  of  Clio  number  of  center  points  k, 
the  central  composite  design  can  be  rude  either  orthogonal  or  uniform  pre- 
cision. Box  and  Hunter  [7]  give  the  characteristics  of  these  designs  for 
various  sizes  of  it. 

Another  important  class  of  second-order  response  surface  designs  is 
the  equiradinl  designs.  That  is,  N >_  (n+l)(n+2)/2  design  points  are  placed 
at  points  on  an  n-dimensional  hypersphere . Figure  9 illustrates  two  equi- 
radial  designs  for  two  variables.  It  is  important  to  note  that  each  of  these 
designs  is  also  equiangular,  in  that  the  H points  are  placed  in  such  a way 
as  to  form  N equal  angles  from  the  design  center.  Biles  and  Swain  [5]  have 
demonstrated  the  construction  of  such  designs  based  on  orthogonal  n-simplex 
designs,  thus  achieving  an  efficient  rotatable  design. 

Having  estimated  the  m second-order  regression  equations  (42)  and  formu- 
lated the  appropriate  optimization  problem,  it  remains  to  apply  mathematical 
programming  to  obtain  a solution.  For  the  constrained  formulation,  any  of 
the  following  procedures  could  be  employed:  (a)  Box's  complex  search  [8]; 

(b)  Rosen's  gradient  projection  method  [23];  or  (c)  one  of  Zoutendijk's 
methods  of  feasible  directions  [28].  For  the  weighted  objective  function 
formulation,  these  same  three  procedures  are  applicable.  For  the  goal 
programming  formulation,  Ignizio's  [17]  procedure  based  on  the  method  of 
Griffith  and  Stewart  [13]  is  computationally  efficient. 

Figure  10  illustrates  a complex  search  applied  computationally  to  second- 
order  response  surfaces.  Figure  11  gives  a flow  chart  for  a simulation/ 
optimization  method  which  employs  the  computational  version  of  Box's  complex 
search  applied  to  second-order  regression  equations  estimated  from  data 
obtained  from  simulation  trials  performed  at  each  point  in  a second-order 
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response  surface  design. 


i 

I 

i 


response'  surface-  design.  Appendix  C gives  the  program  listing  for  this 
procedure . 

Once  an  "optimal"  solution  has  been  obtained,  it  is  necessary  to  compute-, 
a con£ide>nce  interval,  say  at  90%,  about  the  predicted  "optimum."  That  is, 
we  say  that 

(43) 


p l y -j c 1 'fj  i>’jr'l=  °-9> 


JU 

where  y and  y.  are  lower  and  upper  bounds,  resnectivelv,  on  the  90% 
jn 

confidence  interval  for  response  q”  at  X”.  Draper  and  Smith  [10]  give  the 

procedure  for  computing  this  confidence  interval.  If  the  range  (y  -y . ) 

Ju  .1 

is  found  to  be  excessively  large,  it  may  be  necessary  to  perform  confirmation 
simulations  in  the  vicinity  of  the  predicted  "optimum,"  or  even  an  entire 
repetition  of  the  second-order  design  in  a reduced  experimental  region  about 
the  predicted  "optimum." 


Second-Order  Response  Surface  Anal ys i s o f_ box's  Comp! ex  Search  Experiment s 

Another  approach  to  the  simulation/optimization  problem  is  to  combine 
second-order  response  surface  analysis  with  Box's  complex  search.  Such  a 
procedure  would  he  applicable  to  any  of  the  problem  formulations  described 
earlier.  The  general  procedure  for  this  method  is  as  follows: 

1.  Perform  a set  of  Box's  complex  search  experiment s as  described 
in  the  first  procedure. 

2.  Using  the  simulation  responses  at  each  search  point,  fit 
second-order  regression  equations. 

3.  Apply  Box’s  complex  search  computationally  to  these  equations. 

The  advantages  of  this  two-phase  procedure  over  either  of  its  "component" 
techniques  are  that  (1)  it  provides  an  objective  stopping  mechanism  for  the 
experimental  P.ox  complex  search,  and  (2)  it  extracts  latent  information 
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about  the  unknown  response  functions  r^(X),  j 1 , . . . ,m  that  the  searcti 
alone  would  not  discover.  The  following  is  a description  of  this  simulation/ 
optimization  procedure . 

The  Box's  "complex"  search  phase  of  simulation  experimentation  proceeds 


as  follows: 

1.  Randomly  generate  a set  of  N > n + 2 search  points 
satisfying  the  known  bounds  5.  xj_  ;Lrj>  i = l>>-**n- 

2.  Perform  a simulation  trial  at  each  of  these  N search 

r 

points  X ,k=l,...,N  and  record  the  mN  simulated  responses 


y j , .1=1,..., m,  k=l , 


,N. 


3.  Where  a given  search  point  X is  found  to  violate  a constraint 
of  the  form 


H.  /vK\ 

h niK) 


d . , j=l, . . . ,m 


(6) 


generate  a replacement  point  X , perform  a simulation  trial  at 
XK  , and  assess  the  feasibility  of  these  new  responses  yV  , 
j=l,...,m.  Repeat  this  step  until  exact  X search  points  satisfy- 
ing both  the  known  bounds  and  any  constrained  responses  (6)  art 

I 

obtained. 

A.  Identify  as  x'v  that  search  point  yielding  the  least  desirable 

set  of  responses  y') , j=l,...,m.  Compute  the  centroid  XC  of  the 

w c 

N=1  remaining  search  points.  Project  X through  the  centroid  X 

to  a new  point.  If  this  point  satisfies  the  known  bounds,  per- 

w' 

form  a simulation  trial  to  obtain  the  response  v , j=l,..,,in  at 

Xv;  . If  any  responses  fail  to  satisfy  the  constrained  responses 

(6),  shorten  the  stop  from  Xw  through  X'  , and  repeat  this  step 

w» 

until  a feasible  new  point:  X is  obtained.  This  point  replaces 
XV<  in  the  "complex". 
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5.  Repeat  step  4 until  either  of  the  following  conditions  occurs: 

(a)  an  estimated  solution  is  obtained  at  (X,Y);  or 

(b)  a predetermined  number  of  search  points  K (both  those 
satisfying  (6)  and  those  not),  is  obtained  where 

K.  >_  (n+1)  (n+2)/2. 

When  the  complex  search  phase  of  simulation  experimentation  has  been 
terminated,  a best  available,  feasible  solution  (X,Y)  is  retained  for  future 
reference.  An  examination  is  made  of  the  arrangement  of  search  points 
x\  k=1,...,K  in  the  known  experimental  region  <_  £ c^,  i=l,,.,,n  to 

detect  any  regions  that  have  a sparse  density  of  search  points.  Any  such 
region  will  have  a search  point  placed  at  or  near  its  center  and  a simulation 
trial  performed  there.  This  process  guarantees  coverage  of  the  entire  ex- 
perimental region.  It  is  also  worthwhile  to  compute  the  centroid  of  the  T 
total  search  points  and  to  perform  R replicates  of  a simulation  trial  (with 
different  sets  of  initial  random  number  seeds)  at  this  centroid  point.  This 
total  set  of  W=T+R  simulation  trials  comprises  a random  experimental  design 
in  the  factor  space  n.  < x.  <_  c^,  i=l,...,n  with  something  approaching  the 
minimum  number  of  p-ints  for  estimating  second-order  response  surfaces. 

This  minimum  number  of  points  is  (n+1) (n+2) /2,  corresponding  to  the  number 
of  regression  coefficients  in  the  model  (42).  Multiple  linear  regression 
is  then  employed  to  fit.  each  of  the  m second-order  surfaces  to  the  W sets 
of  responses  y}.  , j=l,...,m,  k=l,...,W.  A Box's  complex  search  procedure  is 
applied  computationally  to  this  system  of  m second-order  response  surfaces 
according  to  the  particular  problem  formulation.  The  solution  obtained  by 
this  procedure  (X',Y")  is  compared  to  that  obtained  by  complex  search  (X,Y) 

If  (X* , Y*)  is  preferred  to  (X,Y)  and  XvV  and  X differ  significantly,  one  or 
more  confirmation  simulation  trials  must  bo  performed  at  (X*,Y*)  and  possibly 
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at  (X,Y)  as  well  to  enable  a choice  of  the  "optimum"  solution. 

Figure  12  illustrates  a second-order  response  surface  fitted  to  a set 
of  points  obtained  via  an  experimental  complex  search,  and  then  a computational 
complex  search  applied  to  this  surface.  Figure  11  gives  a flow  chart  for 
this  simulation/optimization  procedure.  Appendix  D is  a program  listing  of 
this  procedure. 

SIMULATION  TFST  MODLLS 

Three  simulation  test  models  were  used  to  evaluate  each  of  the  four 
optimization  methods.  These  simulation  models  are  as  follows: 

1.  A stochastic  inventory  system  as  described  in  Ignall  [16]  and 
Hunter  and  Naylor  [13]. 

2.  A tank  duel  model  as  described  in  Montgomery  and  Bettencourt  [21], 

3.  A naval  minefield  evaluation  model  as  described  by  Bailey  and 
Weedon  [1]. 

Exploratory  work  was  performed  but  not  completed  on  a four  .h  simulation 
model,  the  SPEARS  anti-aircraft  Naval  defense  model  by  Kaplan  et  al  [IS]. 

These  models  are  discussed  in  detail  in  the  following  sections,  but  test 
results  are  discussed  later  In  tills  report. 

Stochast  ic  Inventory  Syst  em 

Ignall  [161  and  Naylor  and  Hunter  [15]  describe  the  application  of 
experimental  design  to  the  optimization  of  computer  simulation  responses. 

They  employed  as  a simulation  test  model  a discrete-event  simulation  of  a 
stochastic  inventory  system  in  which  mean  daily  demand  and  order  lead  time 
are  random  variables  with  known  probability  density  (mass)  functions. 

The  lone  simulation  response  was 

y = mean  daily  cost,  $ 
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which  was  tin.1  sum  of  carrying,  setup  and  shortage  costs.  Tlio  two 
controllable  variables  wore 

x^  = reorder  point  (ROP)  , 

x^  - economic  order  quantity  (KOQ), 

The  optimization  problem  was 


minimize  y = n(xj  ,x,() 


subject  to 


-5  <_  x^  <_  90 


5°  £ \2  < 250 


Ignall  [16]  found  as  a solution  y = $76  at  Xj  = 45  and  X2  = 175  units,  for 
a given  set  of  values  for  the  several  constants  in  the  model. 

The  manner  in  which  this  simulation  test  model  was  utilized  for  the 
current  research  did  not  involve  an  actual  simulation  program.  Rather,  the 
20  x 9 response  table  was  used  as  the  simulator.  This  response  table  is  shown 
in  Figure  13.  For  a given  (x^^),  an  interpolation  was  performed  in  this 
response  table  to  find  v.  This  enabled  a comparison  of  the  four  optimization 
procedures  without  introducing  a confounding  influence  from  the  simulation 
model.  It  should  also  bo  noted  that,  since  the  "simulation"  produced  a 
single  response  v,  this  simulation  test  model  did  not  fully  test  the 
optimization  procedure  but  principally  evaluated  the  interface  between  the 
simulator  and  optimization  modules. 


Tank  Duel  Model 


In  describing  the  application  of  multiple  response  surface  methods 
In  computer  s i no lat  ion , Montgomery  and  Bettencourt  [21]  employed  a stochastic 
simulation  od.l  of  a tank  duel.  The  mol  simulates  brief  fire  engagements 


between  two  armored  vohicl 


A stationary  defending,  vehicle  (Blue  Tank) 


I 

i 

I 
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fires  first  at  a f u 1 ly-exposcd  attacking  vehicle  (Red  Tank).  The  engage- 
ment ends  when  a kill  occurs  or  a predetermined  time  limit  of  120  seconds 
expires . 


The  input  variables  to  the  tank  duel  model  are  presented  in  detail 

in  [21].  But  for  purposes  of  evaluating  and  comparing  the  four  simulation/ 

optimization  procedures  described  in  this  report,  the  following  independent 

varicibles  x.  and  response  variables  y.  were  chosen: 
i '3 


X2 

yJ. 


mean  time  to  fire  first  round  for  the  Blue  crew  (sec.), 
mean  time  between  rounds  (see.), 
probability  of  Blue  victory, 


expected  number  of  rounds  fired  by  the  Blue  tank. 


The  optimization  problem  was  framed  as  one  of  constrained 
optimization,  as  follows: 

maximize  y = g^(x^,x.^) 

subject  to  5 £ x^  <_  25  5 £ x2  £ 25  0 _<  y2  2 


Mine  Hunt  in  ; [ 1 ] 

This  model,  written  in  FORTRAN  IV,  is  a Monte  Carlo  computer  simulation 
which  is  used  to  evaluate  mine  hunting  tactics.  This  program  may  be  used 
either  to  evaluate  minefields  or,  in  conjunction  with  a sustained  attrition 
minefield  evaluation  model  , to  compare  the  effectiveness  of  minesweeping  to 
that  of  mine  hunting.  The  basic  operation  is  to  move  a mine  hunting  ship 
through  a minefield.  Each  Lime  the  ship  passes  a mine  or  minelike  object, 
its  distance  from  that  object  is  calculated,  and  a random  number  in  the 
(0,1)  range  is  drawn  and  compared  with  the.  probability  of  the  ship's  sonar 
detecting  the  object.  If  the  object  is  detected,  a similar  procedure  determines 


if  I he  object  i ••  cor/ectly  classified  and,  in  the  case  of  mines,  neutralized. 
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The  ship's  course  is  altered  to  avoid  mines  it  lias  detected.  Actuation  and 
detonation  of  mines  and  damage  of  ships  is  simulated  for  both  mine  hunters 
and  traffic  ships. 

Inputs  include  the  number  of  mines,  shape  of  minefields,  density  of 
bottom  clutter,  number  of  hunters,  length  of  operation,  effectiveness  of  sonar, 
condition  of  bottom,  number  of  traffic  ships,  and  method  of  hunting.  Geo- 
graphically, the  model  is  limited  to  an  area  20  miles  square  (a  function  of  the 
number  of  mines)  while  the  force  limits  are  100  hunters,  5 traffic  ship 
types,  60  types  of  mines  or  minelike  objects  (e.g.,  bottom  clutter),  and 
750  individual  mines  or  minelike  objects. 

Outputs  include  mines  detected,  mines  neutralized,  ship  damage,  and 
threat  of  tin  minefield  as  a function  of  time. 

Two  separate  problems  were  used  to  evaluate  the  four  optimization 
procedures  with  the  miiu  hunting  model.  In  each  of  these  problems,  the  op- 
timization was  sought  from  the  standpoint  of  the  force  which  deploys  the 
minefield.  That  is,  the  objective  was  to  maximize  minefield  effectiveness 
in  terms  of  damage  to  ships  which  attempt  to  traverse  the  field. 

Problem  1 used  two  mine  types,  arming  delays  and  a total  of  58  mines. 
Problem  2 used  only  one  mine  typo,  no  arming  delays  and  oi  mines.  In  Problem 
1,  all  traffic  ships  (i.c.,  the  target  ships  as  opposed  to  the  mine  sweepers) 
attempted  to  transit  the  mine  field  during  the  period  119-120  hours,  while 
in  Problem  2 traffic  was  uniform  over  the  entire  120  hour  period. 

Additional  information  about  the  problem  that  was  common  to  both  plans 
is  as  follows: 

• The  minefield  was  to  cover  an  area  one  nautical  mile  wide  and  ten 
nautical  miles  long.  (The  input  provided  was  for  a .1800  ft.  x .0  nm. 
field;  this  area  required  fewer  mines,  lienee  reduced  running  lime. 
It  was  assumed  that  mine  requirements  were  proportional  to  the  area 
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of  tlie  minefield  so  that  the  pj  an  for  a Innt.  x .10  nra.  area  used 
(6076.1/1S00)  x the  number  of  mines  in  the  1800'  x 10  nm.  area). 

• Water  depth  was  constant  (hence  for  a given  mine  type/target  type/ 
target  speed  combination,  the  characteristics  of  a mine-target 
interaction  were  the  same  everywhere) . 

• The  target  ships  travelled  at  three  knots  with  a navigational 
error  of  100  yards  (i.e,,  the  perpendicular  distance  from  the 
target's  actual  location  to  the  intended  location  was  normally 
distributed  with  a mean  of  zero  and  a standard  deviation  of  100 
yards) . 

• Two  minesweepers  were  available  and  had  a combined  capability  of 
making  ten  transits  per  day  at  a safe  speed  of  five  knots.  They 
were  assumed  to  have  actuation  characteristics  identical  to  the 
target  but  were  safe  from  damage. 

• The  mines  were  assumed  to  be  uniformly  randomly  distributed 
throughout  the  minefield. 

The  objective  of  each  plan  was  to  provide  a minefield  with  a 90 
percent  chance  of  damaging  (mission  about  level)  at  least  two  of  the  ten 
target  ships. 

The  simulation/optimization  problem  for  each  of  these  two  minefield 
evaluation  problems  was  formulated  ns  follows: 

x^  = mean  of  the  arming  time  distribution,  where  arming  time  was 
Poisson  distributed;  6 £ x^  < 12  lirs. 

X£  ~ mean  of  the  mine  count  distribution,  where  mine  count  was 
Poisson  distributed;  3 < x^  < 5 counts. 

y^  = average  ship  damage,  ships  (summed  over  all  ship  typos) 


i « 7 


, 6 average  "threat"  in  period  j-1,  where  period  1 is 


yj  ’ • 

0 - 24  hrs.,  period  2 is  24-48  lira.,  etc. 


42 


A constrained  optimiznt ion  problem  was 

max  y1 

subject  to  6 £ x^  £ 12 ; 3 < x.^  £5;  0 _<  y ^ £0.25,  j = 2,  ...»  6. 

The  simulation/optimization  interface  for  the  minefield  evaluation 
model  involved  the  following: 

• The  main  program  was  modified  to  become  Subroutine  S1MUL  in  the 
optimization  program. 

• In  Subroutine  S1MFL.  each  mine  was  assigned  random  arming  tine 
and  a random  ship  count  from  Poisson  distributions  having  para- 
meters x^  and  x,,  respectively. 

• In  Subroutine  AVGTHR,  in  which  the  damage  to  traffic  ships  is 
tallied  and  minefield  "threat"  is  computed,  special  program 
segments  were  written  to  give  values  to  the  response  variables 
y y j = 1 . • • • , 6 • 

Air  Attack  on  a Surface  1*1  eet  [18] 

This  model,  written  in  FORTRAN  IV,  evaluates  antiair  warfare  tactics 
related  to  the  following  factors: 

i.  Force  disposition  - e.g.,  dispersed  or  integral  formation  concepts 

ii.  Fire  control  doctrine  - e.g.,  shoot,  look,  shoot; 

iii.  Command  and  control  - e.g.,  coordinated  or  autonomous  assignment 
doctrine ; 

iv.  Cover  and  deception  - e.g., "look  alike"  ships  and  decoy  forces; 

v.  Deployment  of  remote  sensors  - e.g.,  ant ielc ctronic  warfare  air- 
craft and  radar  picket  ships; 

The  effects  of  enemy  force  factlcs  evaluated  include  multidirectional, 
multialt i t ud  Inal  at  tacks,  use  of  electronic  countermeasures  and  penetration 
aids,  "rol.1  back'  attack  tactics,  and  weapon  loading  and  launch  criteria. 


The  node!  May  be  described  as  a one-sided,  event-store,  Monte  Carlo 
computer  simulation.  Enemy  attack  tactics  are  fully  defined  at  the  outset  of 
each  game  situation  by  inputs  to  the  raid  generator,  and  the  attack  is  conducted 
as  defined  regardless  of  the  evolution  of  the  battle.  Thus,  enemy  attack 
tactics  fall  into  the  category  of  the  "uncontrollable  factors"  k = 1,  ...,  p 
in  equation  (1).  The  defense  force,  typically  a carrier  task  force,  includes 
a surveillance  radar  network  (shipborne  and  airborne  air  search  radars,  both 
conventional  periodic  scan  types,  and  multimode  aray  radars),  a communications 
network  for  track  passing  and  firepower  coordination,  and  defensive  weapon 
systems  (primarily  SAM  batteries  and  their  associated  tracking  radars  and 
fire-control  computers).  Task  force  runs  may  be  made  to  establish  sensitivity 
to  changes  in  ship  configurations  (alternate  radars  or  weapon  systems)  or 
in  ship  dispositions  and  assignment  doctrine.  The  model  was  designed  to  be,  to 
t ht'  greatest  extent  practicable,  machine  and  installation  independent. 

Input  data  to  this  model  define  the  size,  composition,  disposition,  and 
fire-control  tactics  for  both  the  offense  and  the  defense  units.  Enemy  forces 
arc  most  often  composed  of  aircraft  and  submarines  that  launch  antiship  cruise 
missiles  from  stand-off  ranges,  although  other  forms  of  attack  arc  possible 
(e.g,,  gravity  bombs).  The  model  has  the  capability  of  handling  up  to  31 
defense  units  (ships  and  AEW  aircraft)  and  up  to  253  offense  units  (launch 
vehicles,  stand-off  jammers,  cruise  missiles,  etc.).  Numerous  other  limits 
define  the  composition  of  each  major  unit  (e.g.,  radars  pership,  missile  launchers 
per  ship,  magazine  capacity  per  launcher,  cruise  missiles  per  launch  platform, 
flight  path  legs  per  target).  Certain  input  parameters  select  tactical  options 
(e.g..  si n gl e or  dual  missile  salvos)  and  simulation  algorithm  options  (e.g., 
simple  or  complex  radar  detection  models). 

A summary  of  the  processed  input  values  is  printed  at  the  outset  of 
a game  sit  uni  Ion.  This  may  opt  ionally  be  followed  by  a detailed  listing  of 
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DISCUSSION  or  RESULTS 

This  section  discusses  the  results  obtained  from  simulation/optimization 
runs  with  the  four  optimization  procedures  in  combination  with  each  of  the 
three  simulation  models.  The  discussion  is  presented  in  three  segments,  one 
for  each  simulation  model. 

Stochast ic  1 nventorv  _Modc  l [ 15,16] 

Table  1 presents  the  results  of  several  optimization  attempts  with  the 

stochastic  inventory  model.  The  complex  search  procedure  and  the  second- 

order  reponse  surface  approach  of  employing  multiple  regression  to  fit  a 

quadratic  model  to  a set  of  complex  search  points  gave  very  similar  results, 

yielding  near-optimal  solutions  with  about  7n  simulation  trials.  The  second- 

order  response  surface  method  employing  a central  composite  design  produced 

a predicted  optimum  at  a point  for  which  the  actual  solution  was  4 percent 

removed  from  the  true  optimum,  using  only  9 simulation  trials.  The  sequential 

2 

fl-  -order  response  surface  approach,  employing  either  a 2 factorial  design 
o'  i-vortex  simplex  design  (each  augmented  by  a centroid  point),  with  a golden 
section  search  along  the  gradient  direction,  required  about  20n  simulation 
trials  to  yield  near-optimal  solutions.  Among  the  four  optimization  methods 
evaluated  with  this  model,  only  the  first-order  response  surface  method  fared 
poorl y . 

It  is  notable  that,  when  a second-order  response  surface  is  fitted  to 
t lie  14  searcli  points  obtained  from  a Hox  complex  search,  the  predicted  op- 
timum solution  is  actually  slightly  worse  than  the  best  solution  obtained 
via  the  search  itself. 
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Table  1.  Summary  of  Simulation/Optimization  Results 
for  the  Stochastic  Inventory  Model  [15,16] 


Opt imizat ion 

Method 

Starting 

Seed 

Est imated 

X1  X2 

SoJ  ut ion 

y 

Number  of 
Trials 

Complex 

12471 

50 

177 

$76.20 

13 

Search 

21437 

43 

245 

81.00 

14 

First-Order 

Factorial  Design 

17332 

70 

125 

$79.00 

40 

First-Order 

Simplex  Design 

17332 

42 

159 

$76.78 

41 

Second-Order 

Central  Composite  Design 

35188 

53 

183 

$70.32 

$79.00 

(Predicted) 

(Actual) 

9 

Second-Order 

Complex  Search 

14271 

41 

163 

$76.78 

77.52 

(Search) 

(Predicted) 

14 

Notes : 

(1) 

Known  Solution 

X1  = X2  = 175  y 

= $76.00 

(2) 

CPU  times  less 
an  IBM  370/168 

than  5 sec.  per  optimization  run 
computer . 

i on 

(3) 

Complex  search 
were  conducted 

terminated  because  2n = 4 search 
without  an  improved  solution. 

points 

47 

Tank  Duel  Model.  [21] 

Table  2 presents  the  results  of  four  simulation/optimization  runs  with 
the  tank  duel  model.  An  optimum  solution  lies  near  (x^  = 8.2,  X£  = 12.5), 
Producing  a probability  of  Blue  Victory  of  0.61  with  an  expected  firing  of 
2 Blue  rounds.  Thus,  the  "expected  rounds  fired"  response  appears  to  bound 
the  solution.  As  with  the  stochastic  inventory  model,  only  the  sequential 
first-order  response  surface  approach  produces  unacceptably  costly  experimen- 
tation . 

In  contrast  to  the  result  with  the  stochastic  inventory  model,  the  pre- 
dicted solution  obtained  after  fitting  second-order  response  surfaces  to 
the  two  sets  of  responses  obtained  in  the  Box  complex  search  represents  a 
significant  improvement  over  the  best  point  observed  in  the  search. 


Mine  Hunting  Model [ lj_ 

As  seen  in  Table  3,  only  three  optimization  methods  were  evaluated  with 
the  mine  hunting  simulation  model.  The  complex  search  method,  as  with  the 
previous  models,  was  terminated  after  2n  simulation  trials  had  been  performed 
without  obtaining  an  improved  solution.  For  example,  with  problem  1,  the 
estimated  optimum  solution  was  observed  with  search  point  12,  and  since  four 
more  search  points  failed  to  improve  on  that  result  the  search  was  terminated 
after  point  16. 

Focusing  on  the  results  of  the  second-order  response  surface  method,  em- 
ploying a central  composite  design  with  16  simulation  trials,  a solution 
x^  ~ 10.6  hours  mean  arming  delay  per  mine  and  X£  = 3.1  ships  mean  ship 
count  per  mine  results  in  an  expected  kill  of  0.51  ships.  The  threat  pro- 
file is  essentially  zero  during  the  first  72  hours  of  minefield  operation, 
and  Increases  to  levels  between  0.15  and  0.2  during  the  last  48  hours  of 
operation  when  enemy  ship  traffic  actually  attempts  to  traverse  the  field. 


A 
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Tabic  2.  Summary  of  Simula t i on/Opt imizat ion  Results 
for  the  Tank  Duel  Model  [21] 


Optimizat ion 
Method 

Estimated 

X1  X2 

Solution 

yl  y2 

Number  of 
Trials 

Complex  Search 

12.9 

10.1 

0.59 

1.97 

20 

First-Order 

Simplex  Design 

S.2 

12.4 

0.60 

1.97 

55 

Second-Order 
Central  Composite 

8.0 

12.5 

0.61 

2.00 

9 

Second-Order 

13.5 

11.3 

0.57 

1.86  (Search) 

15 

Complex  Search 

8.2 

12.2 

0.61 

1.99  (Predicted) 

Notes:  (1)  Known 

optimum 

at  = 

8.2  see 

. , ?:9  =12.5  sec.  , 

y\  = 0.61, 

y2  = 2 rds . 


(2)  CPU  times  less  than  5 see.  per  optimization  run  on  an  IBM 
370/168  computer. 

(3)  Complex  search  was  terminated  because  2n = 4 search  points 
were  evaluated  without  obtaining  an  improved  solution. 
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This  result  points  out  a difficulty  in  interfacing  an  optimization  procedure 
with  a simulation  model  when  a gaming  strategy  is  involved.  In  a sense,  the 
optimization  procedure  "learns"  the  enemy  tactic  as  the  scries  of  simulation 
runs  progresses.  In  this  instance,  the  optimization  procedure  "sees"  that 
traffic  ships  traverse  the  field  late  in  the  game.  To  counter  this  "learn- 
ing" behavior,  it  is  necessary  to  level  the  threat  over  the  entire  period 
of  minefield  operation.  Since  such  a narrow  feasible  region  would  cause 
complex  search  to  yield  a high  proportion  of  infeasible  responses,  the  pre- 
ferred optimization  method  would  be  a second-order  response  surface  approach, 
either  by  a designed  experiment  or  by  fitting  a quadratic  surface  to  Box's 
complex  search  results. 


SUMMARY  AND  CONCLUSIONS 

This  research  has  investigated  four  alternative  approaches  to  optimization 
of  simulated  systems  having  multiple  independent  variables  x^,  i=l , ...,  n 
and  multiple  simulation  responses  y , j=l , . . . , m.  Of  these  four  methods. 
Box's  complex  search  [8],  a second-order  response  surface  approach  employ- 
ing a central  composite  experimental  design  [6,7]  followed  by  a computation- 
al complex  search,  and  fitting  i second-order  response  surface  to  a suc- 
cession of  complex  search  points,  all  yield  "good"  solutions  in  an  "econom- 
ical" number  of  computer  simulation  trials.  The  sequential  first-order 
response  surface  method  converged  to  near-optional  solutions,  but  required 
an  excessive  number  of  simulation  trials.  This  excessive  number  of  trials 
derived  from  the  golden  section  line  search  along  the  estimated  gradient 
direction.  Had  a polynomial  regression  line  search  employing,  say,  five 
trials  boon  conducted  as  in  Biles  [2,3],  two  complete  cycles  of  gradient- 
determining  blocks  and  line  searches  could  hive  been  performed  in  about  Sn 
to  J On  trials. 
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The  most  crucial  aspects  of  the  experimental  Box  complex  search  ap- 
proach to  siraulat ion/optimization  are  twofold: 

] . A termination  criterion  based  on  either  (a)  a maximum  number  of 
simulation  trials  or  (b)  a certain  number  of  simulation  trials 
without  obtaining  an  improved  solution.  Tiie  typical  criterion 
used  to  terminate  a computational  Box  complex  search  is  con- 


vergence; such  as  Such  a criterion  can  lead 

to  excessive  simulation  trials  in  the  experimental  mode,  however, 


unless  e is  coarse. 


2.  In  an  experimental  region  in  which  the  implicit  constraints  on 

simulation  responses  y. , ...,  y describe  a small  feasible  region, 

l m 

a number  of  infeasible  simulation  trials  will  be  conducted.  In 
fact,  the  probability  that  an  early  simulation  tria.'  is  feasible 
with  respect  to  implicit  constraints  is  P^.  = R(V^.)/R(E),  where 
Pj.  = probability  of  a feasible  trial 
R(Yj.)  = region  subtended  by  the  implicit  constraints 
R(E)  = region  subtended  by  the  explicit  constraints 
Once  n+2  feasible  simulation  trials  have  been  found,  and  complex 
search  begins  to  seek  improved  solutions,  there  is  a higher  prob- 


ability of  obtaining  a feasible  simulation  trial.  As  the  "complex" 
becomes  tighter,  this  probability  approaches  one.  But  the  pos- 
sibility of  a large  number  of  infeasible  trials  early  in  tiie  search 
is  perplexing. 

But  the  capability  of  fitting  quadratic  response  surfaces  to  a set  of 
at  least  (n+l)(n+2)/2  complex  search  points  (this  being  the  number  of  re- 
gress Lon  coefficients  in  the  quadratic  model)  in  an  "incomplete"  complex 

y . 


A 
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The  second-order  response  surface  approach  has  the  advantage  of  re- 
quiring relatively  few  simulation  trials,  but  the  accuracy  of  the  predicted 

A A 

optimum  solution  (X,Y)  depends  on  how  well  the  true  surfaces  are  represented 
with  quadratic  estimators.  In  those  cases  where  there  is  some  knowledge  of 
the  nature  of  the  true  surfaces,  the  modeler  is  better  able  to  hypothesize 
the  fitted  model  (which  may  not  be  quadratic).  Although  the  central  com- 
posite design  has  superior  statistical  properties  to  the  purely  random  de- 
sign, such  as  obtained  from  a set  of  complex  search  points,  the  complex  search 
approach  will  usually  place  more  than  one  simulation  trial  in  the  neighborhood 
of  the  predicted  optimum.  This  affords  a smaller  confidence  region  about  the 
predicted  optimum. 

The  computer  programs  listed  in  Appendices  A through  D offer  a means  of 
interfacing  these  optimization  procedures  with  simulation  models  of  interest. 

Minor  alterations  in  the  optimization  program  would  be  necessary  for  a particular 
simulation  model,  but  major  modifications  in  the  simulation  model  itself  could 
be  entailed.  In  general,  one  could  expect  to  have  to  "custom  fit"  the  optimization 
program  to  a given  simulation  model,  this  is  especially  true  of  the  typical 
naval  simulation  model,  because  modelers  tend  to  "start  from  scratch"  in  dev- 
eloping models  for  their  own  needs. 

This  research  has  provoked  two  recommendations  relative  to  naval  simulation 
modeling  activities  in  the  several  naval  laboratories: 

1.  There  is  a definite  need  to  have  naval  simulationists  apply  up-to-date 
statistical  methodology,  particularly  variance  reduction,  in  their 
modeling  efforts. 

2.  There  is  a need  to  have  naval  simulationists  become  familiar  with 
general  simulation  languages  such  as  GASP-IV  and  SIMSCR1PT  2.5.  It 
is  also  apparent  that  certain  specialized  modeling  techniques,  such 
■is  Q-  cun  , is 

I 


vastly  under  utilized  in  naval  operations  modeling. 


Developing  FORTRAN  simulation  models  from  "scratch"  is  sometimes  necessary, 
but  very  often  considerable  time,  effort  and  expense  is  spent  in  duplicating 
the  kinds  of  functional  capability  contained  in  available  and  tested  simulation 
soltware . 

The  most  pressing  need  in  terms  of  continuing  research  in  the  area  of 
optimization  of  naval  simulation  models  is  to  incorporate  the  techniques  of 
game  theory  into  the  optimization  approach.  So  many  naval  models  embody  a 
"friendly  force/enemy  force"  focus,  with  decision-making  capability  on  both 
sides,  that  optimization  must  entail  game  theory  considerations. 
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APPENDIX  A 


PROGRAM  LISTING 
FOR  AN 

EXPERIMENTAL  BOX  COMPLEX  SEARCH 
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APPENDIX  D 

PROGRAM  LISTING 
FOR  A 

SECOND-ORDER  RESPONSE  SURFACE  METHOD 
BASED  ON  BOX'S  COMPLEX  SEARCH 
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